The phase diagram of water from quantum simulations 
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The phase diagram of water has been calculated for the TIP4PQ/2005 model, an empirical rigid non-polarisable model. The path 
integral Monte Carlo technique was used, permitting the incorporation of nuclear quantum effects. The coexistence lines were 
traced out using the Gibbs-Duhem integration method, once having calculated the free energies of the liquid and solid phases in 
the quantum limit, which were obtained via thermodynamic integration from the classical value by scaling the mass of the water 
molecule. The resulting phase diagram is qualitatively correct, being displaced to lower temperatures by 15-20K. It is found that 
the influence of nuclear quantum effects are correlated to the tetrahedral order parameter. 



1 INTRODUCTION 

Ever since the monumental work undertaken by Bridgman in 
1912^ there has been intense and continued interest in the 
phase diagram of water^. The prediction of the phase di- 
agram serves as a severe test for any model of water 4 £. Al- 
though the first computer simulations of water were performed 
in 1969 by Barker and Watts 6 and 1971 by Rahman and Still- 
inger 7 , the calculation of the complete phase diagram was 
only recently undertaken, using the classical models TIP4P 
and SPC/E 8 . Although the TIP4P model provided a qualita- 
tively correct phase diagram, there was room for improvement 
(i.e. the melting point of ice Ih was situated at around 230K). 
Consequently a new re-parameterisation, named TIP4P/2005, 
was proposed 9 leading to a satisfactory description of a num- 
ber of properties of wate r 10 i u , The TIP4P/2005 model is 
a rigid non-polarisable model designed for classical simula- 
tions. In an indirect fashion TIP4P/2005 implicitly incorpo- 
rates nuclear quantum effects, at least at moderate to high 
temperatures. However, the model fails when it comes to de- 
scribing the equation of state at low temperature si^ or Cp^. 
The origin of the failure is the use of classical simulations to 
describe the properties of water. Quantum effects are present 
in water^ 4 ^— even at "high" temperatures, due to the partic- 
ularly small moment of rotational inertia, engendered by the 
low mass of hydrogen, in conjunction with the relatively high 
strength of the intermolecular hydrogen bonds. 

Nuclear quantum effects can be incorporated into con- 
densed matter simulations via the path integral technique pro- 
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posed by Feynman 20 (for an excellent review see 21 ). Barker 22 
and Chandler and Wolynes 23 showed that the formalism of 
Feynman is equivalent, or "isomorphic", to performing clas- 
sical simulations of a modified system where each molecule 
is replaced by a polymeric ring composed of P beads. The 
TIP4P/2005 model, successful for classical simulation s 24 ! 25 , 
was recently adjusted for use in such quantum simulations 
(the charge located on the hydrogen atom was increased by 
0.02e so as to maintain the same internal energy in a quan- 
tum simulation as the TIP4P/2005 model in a classical simu- 
lation) becoming the TIP4PQ/2005 model 12 . This new variant 
of TIP4P/2005 has been successful in describing the tempera- 
ture of maximum density 26 of water and heat capacities^ 3 -. It 
is for this model that we calculate the phase diagram. 

2 METHODS 

2.1 Path integral Monte Carlo 

The partition function, Qn p t, for a system of TV rigid 
molecules in the NpT ensemble is given (except for an arbi- 
trary pre-factor that renders Qn p t dimensionless) by Qn p t = 
f exp(— PpV)Qnvt dV. In the NVT ensemble Qnvt is given 
by: 
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where P is the number of Trotter slices or "replicas" through 
which nuclear quantum effects are introduced (for P = 1 the 
simulations become classical). Each replica, f, of molecule 
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i interacts with the replicas with the same index t of the re- 
maining particles via the inter-molecular potential U, and in- 
teracts with replicas t — 1 and t + 1 of the same molecule i 
through a harmonic potential (connecting the centre of mass of 
the replicas) whose coupling parameter depends on the mass 
of the molecules (M) and on the temperature (j3 = 1/kgT), 
and through a term (p^u ), named the rotational propagator, 
that incorporates the quantisation of the rotation and which 
depends on the relative orientation of replicas t and t + 1 . Pi- 
oneering work was undertaken by Wallqvist and Berne—, and 
by Rossky and co-workers 28 who used an approximate ex- 
pression for the rotational propagator of an asymmetric top 
(i.e. water). Another technique is that of the stereographic 
projection path integral 29 which has been used to study TIP4P 
clusters^. In 1996 Miiser and Berne 31 provided an expres- 
sion of the rotational propagator for spherical and symmetric 
top molecules. Quite recently the authors have extended the 
expression of Miiser and Berne to the case of an asymmetric 
top^ 2 -, which is the case of water. The propagator is a function 
of the relative Euler angles between two contiguous beads and 
of PT. The internal energy E can be calculated through the 
derivative of the logarithm of Qnvt with respect to j3, and is 
given by the sum of the kinetic and potential energy terms, 

E = ^translational + ^rotational +U = K + U. A more Complete 

account concerning path integral Monte Carlo simulations of 
rigid rotors, and their application to water, can be found in the 
article by Noya et al. 33 . 



2.2 Phase diagram calculation 

The determination of the phase diagram of the quantum sys- 
tem is undertaken in several steps. First the classical phase di- 
agram of the TIP4PQ/2005 model was calculated. To do this, 
for each solid phase a reference thermodynamic state is chosen 
and the free energy of the classical system is determined using 
either the Einstein crystal 34 or the Einstein molecule method- 
ologies 35 . For the fluid phase the free energy of the classical 
system is determined at a reference state by transforming the 
TIP4PQ/2005 model into the Lennard- Jones model, for which 
the free energy is well known 36 . The free energy of the clas- 
sical system under distinct thermodynamic conditions can be 
obtained via thermodynamic integration. This permits one to 
determine an initial coexistence point of the classical system 
for each phase transition by imposing the usual condition of 
equal chemical potential for a given T and p. Gibbs-Duhem 
simulations 37 are then performed to trace out the complete 
phase diagram. The procedure has been described in detail 
in—. At the end of this first step the phase diagram of the 
classical system is known. 

In the second step the chemical potential of the quantum 
system is determined at a reference thermodynamic state, 
again for each phase of interest. It is worth describing this 



procedure in some detail. Let us define the excess quantum 
free energy as the free energy difference between the quantum 
system and its classical counterpart at the same T and p: 
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Thus the free energy of the quantum system, G, can be ob- 
tained if the classical and excess contributions are known. 
The free energy of the classical system was determined in 
the first step, so we shall now focus on the evaluation of 
Qex.Q Q ne d e fi nes a parameter, A, whose purpose is to scale 
the mass of the atoms of the molecule of water such that: 
mo = Affloo = m o,o/h' and m# = Am^.o = m^.o/X', where 
mofi and ninfl are the masses of O and H in the molecule 
of water and where A' = 1/A. Thus one can can slide from 
the quantum limit (for which A' =1) to the classical limit (for 
which A' =0) by simply changing the X' parameter. From the 
relationship G = — kT In £>NpT the derivative of the free energy 
with respect to X' can be calculated 39 , obtaining: 



d(G/NkT) 1 / K \ 
dX> ~T'\NkT/ 



(3) 



This is due to the fact that when the mass of all atoms of the 
molecule are scaled by a factor l/X' the total mass, M, is also 
scaled by a factor 1 /X'. The same is true for the eigenvalues of 
the inertia tensor and thus the energies of the asymmetric top 
appearing in the rotational propagator are also scaled by a fac- 
tor A'. The average of the value in the angled brackets should 
be performed for the value of A' of interest. By using Eq. [3] 
in conjunction with the fact that the total kinetic energy (with 
the translational and rotational contributions) is 3NkT for the 
classical and for the quantum system in the limit of infinitely 
heavy molecules, it can be shown that the chemical potential 
of the quantum system jx can obtained from the expression: 
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To determine the integral of Eq. |4] (i.e. G ex & /NkT) it is suf- 
ficient to perform simulations at decreasing values of A', to 
determine the integrand for each considered value of A' and 
subsequently implement a numerical procedure to estimate the 
value of the integral. It follows from Eq. |4] that the differ- 
ence in chemical potential in the quantum system between two 
phases, Afi = jJ-B — Ha, can be obtained as: 
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This expression states that the difference in chemical potential 
between two phases is simply the value of the difference in 
the classical system plus a correction term that accounts for 
the difference in the quantum excess free energies. Thus, after 
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the second step one knows either the chemical potential of the 
quantum system at a reference state, or similarly the difference 
in chemical potential between two phases, again at a reference 
state. 

The third step in the determination of the phase diagram of 
the quantum system requires the determination of one initial 
coexistence point for each coexistence line. By using ther- 
modynamic integration 38 , the free energy of each phase of 
the quantum system is determined as a function of T and p. 
This provides the location of at least one coexistence point be- 
tween each pair of phases by imposing the condition of iden- 
tical chemical potential, p and T between the two phases. 

The fourth and final step is the tracing out of the complete 
coexistence lines thus yielding the phase diagram. This is 
done by using the Gibbs-Duhem simulations, starting from the 
initial coexistence point determined at the end of the third step. 

2.3 Simulation details 

The expression of the rotational propagator is given in—. 
The propagator is composed of an infinite sum over the en- 
ergy levels of the free asymmetric top rotor. In practice the 
summations are truncated; we adopted the criterion that the 
propagator had converged when the absolute difference be- 
tween the value of the propagator for two consecutive val- 
ues of J (normalised so that p (0,0,0) = 1 for both values 
of J) is less than 10~ 6 per point. A grid of one degree for 
each Euler angle was used, and results for intermediate angles 
were obtained by interpolation. Simulations consisted of 360 
molecules for liquid water, 432 molecules for ices Ih and II, 
324 molecules for ice EH, 504 for ice V and 360 for ice VI. 
The algorithm of Buch et al. 40 was used to obtain a proton 
disordered configuration of ices Ih, III, V and VI simultane- 
ously having zero dipole moment and at the same time sat- 
isfying the Bernal-Fowler rules 4 'i 42 . For ice III we addition- 
ally imposed the condition that the selected proton disordered 
configuration presented an internal energy that lies in the cen- 
tre of the energy distribution shown in Figure 2 of 43 . Direct 
Ih-fluid coexistence simulations used about 1000 molecules 
for the classical system and about 600 for the quantum one. 
Free energies of the classical system were obtained using the 
Einstein molecule methodology for a given proton disordered 
configuration and subsequently adding the Pauling 42 entropy 
contribution (— RT1n(3/2)), The methodology has been de- 
scribed in detail elsewhere 38 . The Lennard- Jones part of the 
potential was truncated at 8.5A and long ranged corrections 
were added. Coulombic interactions were treated using Ewald 
sums. For the solid phase anisotropic NpT Monte Carlo simu- 
lations were performed in which each of the sides of the simu- 
lation box were allowed to fluctuate independently. The num- 
ber of replicas, P in the path integral simulations was selected 
for each temperature by imposing that PT be approximately 



1900 ± 300. This choice guarantees that, for a rigid model of 
water, the thermodynamic properties are within two per cent 
of the value obtained as P tends to infinity. In general the sim- 
ulations of this work consisted of 200,000 Monte Carlo cycles, 
where a cycle consists of a trial move per particle (the num- 
ber of particles is equal to NP where N is the number of wa- 
ter molecules) plus a trial volume change in the case of NpT 
simulations. To increase the accuracy, in the determination of 
the excess quantum free energies, four independent runs were 
performed for each value of A'. In Eq. [5] A/x classlcal is zero 
if evaluated at the coexistence T and p of the classical sys- 
tem. Direct coexistence simulations were significantly longer 
(up to 10 million cycles) and Gibbs Duhem simulations were 
typically ten times shorter. Further details of the path integral 
Monte Carlo simulations, for example, the rotational propaga- 
tor, the acceptance criteria within the Markov chain, the eval- 
uation of relative Euler angles between contiguous beads, can 
be found in 22^. 

3 RESULTS AND DISCUSSION 

To illustrate the methodology used in this work to determine 
the entire phase diagram of the quantum system, we shall de- 
scribe in detail the procedure used to determine the fluid-I^ 
coexistence curve: 

• Determine the melting temperature 7^ lasslcal of the classi- 
cal system at normal pressure. 

• Determine the integral in Eq. H]by performing path in- 
tegral NpT simulations for various values of X' for ice 
Ih and for liquid water. For H2O the values of V = 
1,6/7,5/7,4/7,3/7,2/7,1/7 were used to evaluate this 
integral. 

• Perform thermodynamic integration and determine the 
melting temperature T m of the quantum system at which 
ice Ih and water have the same jj. at normal p. 

• Perform Gibbs-Duhem integration using path integral 
simulations to determine the full Ih-water coexistence 
line. 

Free energy calculations for TIP4PQ/2005 yielded r,£ lassical = 
282K for Ih (p = 1 bar). The same result 282±3K was ob- 
tained from direct coexistence simulations. In direct coexis- 
tence runs 38 ' 44 half of the simulation box is filled with ice and 
the other half with the liquid. NpT simulations at normal pres- 
sure are performed for several temperatures. For temperatures 
above the melting point the ice within the system melts, and 
for temperatures below the melting point the ice phase is seen 
to grow. 

We then proceeded to calculate the integral in Eq. |5]at nor- 
mal p and T — 282K (where the two phases have the same 
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chemical potential in the classical system). To do this the mass 
of the atoms in the TIP4PQ/2005 "molecule" were incremen- 
tally increased by a factor of A. Such a scaling only modifies 
the total mass of the molecule M and the eigenvalues of the in- 
ertia tensor (i.e. the principal moments of inertia), but leaves 
the geometry of the model and the location of the centre of 
mass unchanged. Seven scaling factors between A = 1 and 
A = 7 were used, and the corresponding exact rotational prop- 
agator at T = 282K was calculated for each. Beyond A = 7 
the calculation of the propagator becomes prohibitively ex- 
pensive to calculate, and a large number of simulations would 
be required to reduce the errors to an acceptable level. Path 
integral simulations were then performed, and the integrand 
of Eq. [5] is determined (see Figure la). The integral over the 
curve formed by these results provides the difference in excess 
quantum Gibbs energy between ice I], and water, and therefore 
the difference in free energies between the two phases in the 
quantum system (for T = 282K and p — 1 bar the two phases 
have the same chemical potential in the classical system). As 
can be seen in Figure la the integrand for the liquid-It, calcula- 
tion is reasonably smooth and forms an almost horizontal line. 
For this reason it seems reasonable to extrapolate the integrand 
for A' < 1/7 from the values obtained for larger A'. 

It can be seen that the kinetic energy of ice Ih is higher 
than that of water at 282K and lbar, indicating that nuclear 
quantum effects are significantly larger in the ice Ih phase 45 . 
For rigid models nuclear quantum effects are related to the 
strength of the intermolecular interactions which, in the case 
of water, is dominated by hydrogen bonds. In ice Ih each 
molecule forms four hydrogen bonds with the first nearest 
neighbours, whereas in the liquid phase this number is some- 
what smaller. The more "localised" character of the molecular 
libration in ice Ih with respect to the liquid leads to the higher 
kinetic energies observed. From the results shown in Figure 
la it follows that jj, of ice Ih in the quantum system is 0.1 8&T 
higher than that of water at 282K and lbar, indicating that the 
melting point of ice Ih in the quantum system is lower than that 
of the classical system. By using thermodynamic integration 
at constant p we found that at 258K the chemical potential of 
the solid and fluid phases become identical. This value is the 
T m of the quantum system. In order to corroborate this result, 
direct coexistence runs of the quantum system were under- 
taken. For the runs performed at T — 266K and T = 262K the 
total energy increases with time and reaches a plateau indicat- 
ing the complete melting of the ice slab. At T = 240K we saw 
a slow growth of the ice phase (See|2). At T = 252K the en- 
ergy was approximately constant along the run. This indicates 
that the melting point lies between T = 252K and T = 262K, 
thus we shall adopt the intermediate value of 257K (±5K) in 
agreement with the free energy result of 258K. 

The entire Ih-water coexistence line is obtained by using 
the Gibbs-Duhem integration method. The Gibbs-Duhem in- 



tegration method consists of a numerical integration of the 
Clapeyron equation and requires simply the knowledge of 
the enthalpy and volume difference between the two coex- 
isting phases. The enthalpy of each phase is obtained from 
H = K + U + pV. Note that for each new temperature in the 
Gibbs-Duhem integration a propagator matrix must be calcu- 
lated as the propagator depends on the value of PT, It can be 
seen that the Ih-water curve (0 is essentially parallel to the 
experimental curve, shifted by ss 15K due to the lower melt- 
ing point of the TIP4PQ/2005 model. The aforementioned 
methodology for Ih-water was applied to the remaining phase 
equilibria, leading to the complete phase diagram. A plot of 
the integrand of Eq. [5]with respect to A ' is shown in Figure 1 a. 
The integral of these functions leads to the difference in ex- 
cess quantum free energy between the two considered phases 
in units of NkT. As can be seen the integrand is rather smooth, 
and in most of the cases it can be well described by a straight 
line. It is worth noting that if this were always the case then 
one could obtain a reasonable estimate of the integral simply 
by obtaining the value of AK/ (X'NkT) at A' = 1 /2. In the 
phase diagram of water as obtained from quantum simulations 
of the TIP4PQ/2005 model of water is presented and com- 
pared to the experimental phase diagram. One can see that the 
diagram is qualitatively correct, each phase is situated in the 
correct relation to the other phases. Furthermore, the gradients 
of the coexistence curves are also acceptable in comparison 
to experimental results. The most notable discrepancy is an 
overall shift of 15-20K in the diagram to lower temperatures. 
In |4] the changes in volume along phase transitions obtained 
from the simulations are presented. It can be seen that they 
compare favourably with the experimental results obtained by 
Bridgman in 1912 1 Making use of a recent publication by 
Loerting et al. 46 where the densities of ices Ih, II, III and V in 
the range 77-87K at normal pressure were determined using a 
methodology known as cryoflotation, we decided to study the 
densities at the intermediate temperature of 82K. Additionally 
we also considered a couple of states for the liquid, and one 
for ice VI at room temperature. The results are summarised in 
Q] In general the simulation and experimental results coincide 
to within 1%. Classical simulations of TIP4P/2005 tend to 
overestimate the experimental densities (at 82K) by more than 
3%—. Thus a quantum treatment is absolutely essential if one 
wishes to describe experimental results at low temperatures. 
The results of this table can also be used to estimate the value 
of the transition pressures at zero Kelvin as shown by Whal- 
ley^Z. The estimates obtained in this way were consistent 
with those obtained from extrapolations to zero Kelvin of the 
coexistence lines. The maximum deviation found between the 
two methodologies was «700 bar, which is reasonable taking 
into account the combined uncertainty of all the calculations. 

In order to highlight the differences between quantum and 
classical results for the phase diagram, the quantum and clas- 
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sical phase diagram of the TIP4PQ/2005 model are superim- 
posed in [5] Although the diagrams are qualitatively similar 
there are certain features of interest that can be observed. In 
the classical phase diagram the melting lines are shifted to 
higher temperatures, and solid-solid transitions are shifted to 
higher pressures (for a given temperature) with respect to the 
quantum phase diagram. Another important difference be- 
tween the classical and quantum phase diagrams is that the 
region of the phase diagram occupied by ice II is significantly 
reduced in the classical treatment. In fact in the classical sys- 
tem the ice II-III transition is shifted to much lower tempera- 
tures and ice II is stable only for temperatures below 80K. This 
shrinking of the ice II phase is consistent with recent findings 
by Habershon and Manolopoulos 48 who found that in classi- 
cal simulations of the q-TIP4P/F model 49 ice III occupies the 
region of stability of ice II. This indicates that nuclear quan- 
tum effects play a significant role in determining the region of 
stability of ice II in the phase diagram of water. 

It would be useful to have a rational guide to understand 
the changes in the phase diagram observed when including 
nuclear quantum effects. For this purpose the integrand of Eq. 
|U which facilitates the determination of the quantum excess 
free energy, is shown in Figure lb at a temperature of 200K. 
The average kinetic energy of a harmonic oscillator of mass 
M = Mo/A' and frequency v = VX'Vq is given by—: 



(K) = 
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Upon performing a Taylor series expansion about X' = one 
obtains: 
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For the rigid water model used in this work, one can describe 
the solid phases by a set of 6N oscillators (i.e. phonons). By 
assuming a unique frequency, as in an Einstein like model, one 
arrives at: 
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(8) 

Thus for the Einstein model the integrand of Eq. 4 is well be- 
haved and has both a finite value and a finite negative slope at 
X' = . The results presented in Figure lb are indeed con- 
sistent with this predicted behaviour. From Figure lb is is 
possible to estimate Vo from either the slope or the intercept, 
obtaining "Einstein" like frequencies between 550 cm" 1 and 
450 cm -1 depending on the phase. These are typical values 
for intermolecular librations in water, which are located be- 
tween 50 cm" 1 and 800 cm" 1 . A recent study has shown that 
one can reproduce the heat capacity of ice Ih using a selection 
of six fundamental frequencies selected from this range—. It 



is worth mentioning that the TIP4PQ/2005 model also does a 
good job of calculating C p when used in path integral simula- 
tions^. The excess quantum free energies are obtained from 
the integration of the results of Figure lb. It is evident that at 
200K nuclear quantum effects significantly influence the free 
energies of the solid phases of water. By comparing the re- 
sults of ice II at both 1 and 4112 bar it is seen that pressure 
increases the magnitude of nuclear quantum effects at a given 
temperature although the increase is small (0.08 in NkT units 
for the considered pressures). The relative ordering in which 
the excess free energy increases is II ~VI, V, III and finally Ih. 

The tetrahedral order parameter, q t 52 ^ , was designed to 
measure the degree of tetrahedral ordering in liquid water: q, 
is defined as: 
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(9) 



where the sum is over the four nearest (oxygen) neighbours 
of the oxygen of the i-th water molecule. The angle 9,- / # is 
the angle formed by the oxygens of molecules j, i and k, oxy- 
gen i forming the vertex of the angle. The tetrahedral order 
parameter has a value of 1 for a perfect tetrahedral network, 
and for an "ideal gas" of oxygen centres. We shall make 
use of this descriptor in order to try to rationalise our results. 
The value of q t of each solid at T=200K was obtained by an- 
nealing the solid structure while keeping the equilibrium unit 
cell of the system. In|6l the excess free energy is plotted as a 
function of q t for the proton disordered ices, namely Ih, III, V 
and VI at a pressure of around 3600 bar and a strong correla- 
tion is evident. Ice II has a large value, q t = 0.83, however, 
the impact of nuclear quantum effects on this ice phase are 
smaller than in the rest of the ices. Since ice II is the only 
proton ordered solid considered in this work, it is clear that 
in this case the fixed relative orientations between molecules 
are playing an important role in determining the magnitude of 
the nuclear quantum effects. It would be useful to evaluate the 
impact of nuclear quantum effects on the fluid phase. At 200K 
the fluid is highly supercooled thus is difficult to evaluate its 
Qex.Q -pjjg diff erence between the total kinetic energy of a 
phase and that of the corresponding classical system under the 
same conditions also provides an estimate of the magnitude of 
nuclear quantum effects. One can see in Q] that at 300K and 
15400 bar the kinetic energy of liquid water is only slightly 
lower than that of ice VI. This indicates that the magnitude of 
nuclear quantum effects in the liquid is smaller than that of 
the ice with smallest nuclear quantum effects, ice VI. This is 
consistent with the low value q, = 0.58 of the tetrahedral order 
parameter found for the fluid phase at 300K and 15400 bar 54 . 
It appears that the importance of nuclear quantum effects in- 
creases as the strength of the intermolecular hydrogen bond- 
ing increases. The strength of the hydrogen bonding seems to 
correlate (with the exception of ice II) with the value of the 



5 



4 CONCLUSION 



tetrahedral order parameter. For example, in ice Ih the first 
four nearest neighbours of a given molecule are located in a 
perfect tetrahedral arrangement which is the optimum situa- 
tion to have a strong hydrogen bond. For the rest of the ices 
(and for water) the four nearest neighbours of a molecule form 
a distorted tetrahedron and therefore the strength of the hy- 
drogen bond should decrease. The greater the strength of the 
hydrogen bond the higher the frequency associated with the li- 
brational mode, and therefore the higher the impact of nuclear 
quantum effects. 

The change of the coexistence pressure for a certain tem- 
perature due to the inclusion of nuclear quantum effects can 
be approximated reasonably well by the expression 



P /^classical 



(10) 



Vb-V a 

where the properties on the right hand side are evaluated at 
^classical- It follows from Eq. [TO] that the impact of nuclear 
quantum effects on a given phase transition depends on the 
difference of the excess quantum free energy between the two 
phases, and on the volume change. The excess free energy 
difference between phases decreases in the following order: 
liquid-It, H-Ih, liquid-Hi, H-HI, liquid- V, III-I h , IH-V, liquid- 
VI and finally II- VI. The impact of nuclear quantum effects 
on a certain phase transition will be small when the volume 
change of the phase transition is large, and large when the 
volume change is small. Volume change along the phase tran- 
sitions of water are presented in |4] Taking these two factors 
into account it is clear that the II-III phase transition is most 
affected by nuclear quantum effects (i.e. the excess free en- 
ergy difference is large and the volume change is small), fol- 
lowed by the melting curves of ices (decreasing in the order 
liquid-Ih, liquid-Ill, liquid- V). This is followed by the transi- 
tions Ih-II, m-V and Ih-m. Finally the liquid- VI and the II- VI 
coexistence lines are those least affected by nuclear quantum 
effects. 

4 CONCLUSION 

This work illustrates that the calculation of the phase diagram 
of water, including nuclear quantum effects, is now feasible, 
although admittedly it is computationally expensive (even for 
the simple model considered in this work 8 CPU's were re- 
quired for about 2 years to obtain the phase diagram pre- 
sented). The impact of nuclear quantum effects on phase tran- 
sitions is significant and can be rationalised in terms of the 
degree of tetrahedral ordering of the different phases and of 
the magnitude of the volume change involved in each phase 
transition. The TIP4PQ/2005 model yields a reasonable pre- 
diction of the experimental phase diagram of water. The sim- 
ulation results are consistent with the Third Law of thermo- 
dynamics and predict rather well the densities of the different 




Fig. 1 (a) Integrand of Eq. l](i.e. (K B - K A )/(X'NkT)) as a 
function of X! for transitions A — B. Key: red line with ♦ is liquid-Ih 
at 282K and p = lbar, magenta line with □ is II- V at 200K and 
p = 41 12bar, blue line with ■ is II- VI at 200K and p = lbar, 
magenta line with A is V-VI at 200K and p = 9505bar, blue line 
with O i s Ih-HI at 200K and p = 3306bar, and the red line with • is 
Ih-II at 200K and p = lbar. Error bars (only shown for liquid-Ij,) 
represent the standard error, (b) Integrand of Eq. 0for several ices at 
200K. Results were obtained using path integral simulations of the 
TIP4PQ/2005 model. The lines correspond to a fit of the simulation 
results to a second order polynomial. Results (from top to bottom) 
correspond to ice \ at 3306bar, ice III at 3306bar, ice V at 41 12bar, 
ice II at 41 12bar, ice II at lbar and ice VI at 1 bar. The integral of 
the curves (from to 1) yields G ex ®/ (NkT) which results in (from 
top to bottom) 3.1 1, 2.91, 2.80, 2.68, 2.60 and 2.59. 



phases over a wide range of temperatures and pressures. With 
some delay with respect to the original contributors, this work 
shows that a simple modification of the water model proposed 
by Bernal and Fowler 41 in 1933, can reproduce reasonably 
well the experimental phase diagram of water determined by 
Bridgmann 1 in 1912, providing results at low temperatures 
consistent with the Third Law first stated by Nernst in 1906. 

In concluding this work it is worth commenting on the use 
of a rigid non-polarisable model to represent water. Naturally 
in reality water is both flexible and polarisable 55 so it goes 
without saying that this work is far from the last word on the 
matter, and the results presented here form only a way-point 
on the long road to obtaining a definitive model of water that 
describes all of the facets of this intriguing molecule. That 
said, path integral simulations of the TIP4PQ/2005 model has 
provided us with the best phase diagram of water calculated 
to-date. 

This work was funded by grants FIS2010-16159 and 
FIS201 0-1 5502 of the Direccion General de Investigacion and 
S2009/ESP-1691-QF-UCM (MODELICO) of the Comunidad 
Autonoma de Madrid. The authors would like to thank Prof. 
J. L. F. Abascal for many helpful discussions. 
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0.010° 2.510 6 5.0 10 6 7.5 10 6 1.0 10 7 

Monte Carlo cycles 

Fig. 2 Plot of the total energy (E) per particle from the liquid-I n 
direct coexistence simulations of the quantum system at p = 1 bar. 




140 185 230 275 320 200 230 260 290 320 

temperature (K) 

Fig. 4 Molar volume change (Av = vb — va) along the phase 
boundaries from (Left) path integral simulations and (Right) 
experimental results 1 . (w indicates the liquid phase). 



12000 




140 160 180 200 220 240 260 280 300 320 
temperature (K) 

Fig. 3 Phase diagram of water from path integral simulations of the 
TIP4PQ/2005 model. Experimental results (blue points) are also 
showni2& 
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Fig. 5 Classical phase diagram of the TIP4PQ/2005 model (dashed 
blue lines) compared to the diagram obtained from path integral 
simulations (solid red lines). 



7 



4 CONCLUSION 
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% 

Fig. 6 Correlation between the excess quantum free energy and the 
tetrahedral order parameter for ices at 200K and p ~ 3600 bar. For 
ice VI the excess quantum free energy at 1 bar was increased by 
0.08NkT to estimate its value at 3600 bar. 
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Table 1 Densities of ices and liquid of water (in g/cm 3 ) under different thermodynamic conditions. All results were obtained for p = lbar, 
except the two labelled with an asterisk for which p = 15400 bar. Experimental densities for ice III correspond to the experimental values for 
ice IX ( fl — , b — ) the proton ordered form of ice III. Experimental results are from Refs. 46 i 46 i 58 'H The average values of U and K obtained 
from simulations (kcal/mol) are also shown. 



Phase 


r(K) 


p (sim.) 


P (exp.) 


U 


K 


Ih 


82 


0.927 


0.932 


-14.302 


1.914 


II 


82 


1.189 


1.211 


-14.136 


1.774 


II 


123 


1.185 


1.190 


-14.046 


1.837 


III 


82 


1.148 


1.169 a ,1.160 fo 


-14.040 


1.863 


V 


82 


1.252 


1.249 


-13.883 


1.808 


VI 


82 


1.335 


1.335 


-13.745 


1.790 


VI(*) 


300 


1.383 


1.391 


-13.055 


2.475 


liquid(*) 


300 


1.312 


1.311 


-11.997 


2.395 


liquid 


300 


0.997 


0.996 


-11.897 


2.366 
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